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Abstract 



The problem of image segmentation is known to become particularly 
challenging in the case of partial occlusion of the object(s) of interest, 
background clutter, and the presence of strong noise. In the case when 
the segmentation is performed by means of active contours, a typical way 
I I to overcome the above difficulties is to impose an a prion model on the 

shape of the contours - a model which can be either analytical or prob- 
abilistic in its nature. Effectively, the model represents some shape con- 
^ straints, which are intended to regularize the segmenting contour in the 

^ case when imagery data alone fails to provide sufficient information for 

I I determination of its optimal configuration. In practice, the shape models 

are typically learned based on training sets of examples of the object(s) 
" ^' of interest. In such cases, the goodness of modeling depend on the size 

J!L of the training set, with more examples leading to more reliable and ac- 

curate models. Unfortunately, when the number of training samples is 

, relatively small, the resulting model can be inadequate, thereby tending 

to bias the segmenting contour towards a suboptimal solution. To over- 
come this deficiency, the present paper introduces a novel approach to 
modeling of shape priors. Specifically, in the proposed method, an ac- 
tive contour is constrained to converge to a configuration at which its 
geometric parameters attain their empirical probability densities closely 
matching the corresponding model densities that are learned based on 
training samples. It is shown through numerical experiments that the 
proposed shape modeling can be regarded as "weak" in the sense that it 
minimally influences the segmentation, which is allowed to be dominated 
by data-related forces. On the other hand, the priors provide sufficient 
constraints to regularize the convergence of segmentation, while requiring 
substantially smaller training sets to yield less biased results as compared 
to the case of PCA-based regularization methods. The main advantages 
of the proposed technique over some existing alternatives is demonstrated 
in a series of experiments. 
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1 Introduction 



Image segmentation is know to be a problem of fundamental importance in nu- 
merous applications of computer vision and image processing yj , some standard 
examples of which include (yet not limited to) medical imaging [2]|3], surveil- 
lance [4j[5], robotics [6], and control |I7 . In all these applications, image seg- 
mentation is employed to partition a data image (or a sequence thereof) into a 
number of its fragments which are associated with different classes, normally an 
object and its background. In 2-D imaging, such a partition can be performed 
using deformable curves (also known as active contours), whose optimal configu- 
ration maximizes a measure of statistical dissimilarity between the segmentation 
classes. Thus, at some fundamental level, the segmentation is based upon the 
fact that, in a properly defined domain of image features, the object and its 
background may have distinct statistical characteristics. Unfortunately, the va- 
lidity of this assumption is compromised in the presence of strong noises and 
partial occlusions, which necessitates the development of robust segmentation 
tools. 

In the past, a variety of classical algorithms [8 ■ 12 were proposed to segment 
objects of interest based on imagery data alone (see also [iSl and references 
therein for a summary of such techniques). Through employing edge-related 
and region-based features of data images, the above methods are capable of 
providing reliable segmentation results on conditions of moderate noises and 
unoccluded objects. Unfortunately, the methods are known to be prone to erro- 
neous segmentation in practical scenarios, in which some parts of the object of 
interest appear to be occluded, missing, or corrupted by strong noises. In these 
situations, a standard remedy is to enhance the image segmentation through 
the introduction of prior shape knowledge. 

A variety of segmentation algorithms incorporating shape priors have been 
proposed in the literature 3p4 ■ 20 . Most of these algorithms poses the segmen- 
tation task as a combination of two competing optimization problems: whilst 
the first problem maximizes the likelihood of contour's configuration based on 
image-related information, the second problem minimizes the extent to which 
the contour violates the shape constraints imposed by the prior model. Thus, 
for example, the prior model in '15] is defined by the average shape of a set 
of training segmentations. Subsequently, the active contour is constrained to 



deviate the least from this average. A similar approach is exploited in 18 , with 



a diiTerent cost function used to assess the above deviation. While useful in the 
situations when the shape of the active contour needs to be strongly constrained 
(as it would be the case with occlusions), the methods of 15 18 can result in 
useless segmentation when the actual shape of the object happens to deviate 
considerably from its sample average. 

In 14 , a parametric shape model is derived by performing principle compo- 



nent analysis (PCA) on a training set of level-set functions 21 . Subsequently 



the level-set function representing the actual active contour is modeled as a lin- 
ear combination of the principal components (also known as eigenshapes 22 ) 
offset by the algebraic mean of the training functions. A slightly different for- 
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mulation of the same idea is presented in [s], which uses a PCA-based shape 
modehng in combination with the region-based active contours of 12 . As will 
be demonstrated later in this paper, one of the main drawbacks of the PCA- 
based shape modeling lies in the dependency of its validness on the size of the 
training sets in use. In particular, small training sets lack the ability to rep- 
resent all possible appearances of the desired object, and hence the resulting 
PCA model could introduce a significant bias in the estimated shape. More- 
over, since it is based on Euclidean- type metrics, PCA could be easily affected 
by misrepresentative training samples (outliers), whose negative influence will 
be particularly destructive in the case of small training sets. Finally, the PCA- 
based methods [3 14 22 require the use of an alignment procedure which must 
be performed at each iteration of the segmentation process. This alignment step 
is posed as an additional minimization problem, which substantially increases 
the overall computational complexity of segmentation. 



In a more recent work reported in 19 , a segmentation method based on 

is 



the concept of distribution tracking has been proposed. The method of 19 
based on learning both the shape and appearance models, and it is designed 
in a way that eliminates the need to compute pixel-wise correspondences for 



alignment purposes. Moreover, 19 alleviates the problem of high dimensionality 
associated with level set functions via performing PCA on a smaller set of control 
points representing the training shapes. Although the above algorithm appears 
to address some of the shortcomings of previous techniques, the reliability of its 
shape model still depends on the size of the training sets in use. Therefore, the 
associated disadvantages remain. 

The method of is based on the earlier results of 23 , in which the concept 
of segmentation via distribution tracking was first described. In particular, the 
idea in [23| is to learn the probability distribution of a photometric parameter 
(or a set thereof) of the object class, followed by propagating the active contour 
towards a configuration at which the image region encompassed by the con- 
tour has the same parameter distributed similarly to the learned model. In this 



regard, the present work extends the idea of 23 to tracking of morphological 
parameters/features. Specifically, given training data, the probability densi- 
ties of a set of morphological and photometric parameters are estimated first. 
Subsequently, the active contour is propagated towards a configuration which 
results in a "match" between both the photometric and morphological distri- 
butions. The latter can be seen as playing the role of a "weak" shape prioi[^ 
which nevertheless is informative enough to provide an effective regularization 
force. In particular, in our experimental study, we will show that the "weak" 
model outperforms the PCA-based shape representation for the case of small 
training sets. Moreover, an additional advantage of the "weak" model stems 
from the fact that the morphological parameter(s) can be chosen to be invari- 
ant under a class of geometrical transformations (such as Euclidean or affine), 
which effectively eliminates the need for intermediate alignment routines. 

^The notion of "weakness" is introduced to stress the fact that many shapes may have an 
identical distribution of their geometric parameters, such as curvature. 
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The rest of the paper is organized as follows. In Section II, a method for im- 
age segmentation via tracking both texture-related and morphological features 
is detailed. Some essential numerical details are provided in Section III. Section 
IV presents a series of experimental results, which demonstrate some principal 
advantages of the proposed technique. Finally, Section V summarizes the paper 
with a discussion, conclusions, and an outline of our future research. 



2 Tracking of distributions 

2.1 Tracking of photometric features 

Let / be a scalar- valued image defined over an open subset e . To make the 
discussion general, the image / can be transformed into a vector-valued image 
J : 17 -> of its associated (local) features, with N equal to the number of 
features. The transition from I to J can be formally described by means of a 
map Ai : 1 J, which is equal to identity in the case when the only features of 
interest are the gray-levels of /. In general, however, at each x S fi, J(x) may 
consist of, for example, the local statistics of /(x), its associated Laws' texture 



features 24 , multi-resolution moments [25 , or any combination thereof. The 
choice of features is wide and diverse; a specific combination of such features is 
usually selected based on the application at hand. 

Let ilm be a subset of fl over which the object of interest is supported. The 
segmentation algorithm proposed in this paper is based on the following two 
assumptions. First, it is assumed that a set of training images is available with 
their corresponding fl„i being identified. Second, it is assumed that for each 
X G flm , the N components of J(x) are realizations of N independent random 
variables. The latter assumption is obviously an approximation, whose validity 



can be ameliorated by means of a whitening transform 26 . 

Let Pm{z), where z e M^, be the joint probability density function (pdf) of 
the image features corresponding to Qm- Due to the assumption of statistical 
independence, Pm(z) can be factorized as Pm(z) = Y\k=i Pmi^k) , with Pm{zk) 
being the pdf of the A:-th image feature Zk € K. Given flm, the pdf Pm{zk) can 
be approximated according to 

J L„ j^(z,-jr(x))rfx | 

pM = £^ IZd^ 1' 

where K denotes a (positive valued, normalized) kernel function, J]^ (x) is the 
value of the A;-th component of a training image J**" at location x, and 8 stands 
for the operator of assemble average taken over the training set. Note that due 
to the normalization of K (i.e. J^K{s)ds = 1), the estimate injl) is, in fact, 
a non-parametric kernel density estimate of the true pdf of [27 . It should 
also be emphasized that the estimates Pm{zk) are assumed to be pre-computed 
before the actual segmentation is initialized. 

Now, let / be an observed image to be segmented, and J be its corresponding 
feature image. Also, let C $7 be an arbitrary (yet non-empty) subset of il. 
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In this case, an approximation similar to that in ([T]) can be employed to estimate 
the empirical pdf P(z | Qm) of the features of / observed over J7i„. Due to the 
assumption of statistical independence, P(z | f2i„) can be factorized similarly 
to Pm, with its N factors computed according to 

L K(zk — Jfc(x)) dx 
p{zk\n,n)^-^ 7. ^ , withA; = l,2,...,iV, (2) 

so that P{z I Qin) = Ylk=lPi^k I ^in)- 

The ultimate goal of image segmentation is to "deform" fi^n so as to make it 
closely approximate the subset $7™ that supports the object of interest. More- 
over, whenever r2j„ coincides with O^, the distributions P{z \ flin) and Pm(z) 
should be close one to another in some sense. In this paper, as a measure of sim- 
ilarity between probability distributions, we use the Bhattacharyya coefficient 



B which is defined as given by 28 ■ 30 



Bin^n) = I y^KMPWhn)dz. (3) 

Alternatively, the statistical independence of Zk allows rewriting the above ex- 
pression in a slightly simplified form 



N N 



B{Vtin) = J]^ Bk{ilin) = n / VPmizk) P{zk \ ^in)dzk. (4) 

7 1 7 1 



fc=l fe=l 



The Bhattacharyya coefficient B{D,in) achieves its maximum value of 1 when 
Pm{zk) ~ p{zk I ^7in),Vfc, which happens when the sets f2i„ and VL„i coincide 
with each other. Conversely, by maximizing i? as a function of r2i„ one can 
reasonably expect to minimize the discrepancy between f2i„ and " the fact 



that forms the basis of the segmentation method of 19 23 . Needless to say, 
the maximization of B over all possible rii„ is a computationally intractable 
problem, which happens to have an elegant solution when rii„ is replaced by its 
implicit definition in terms of a level-set function as explained next. 

2.2 Level-set formulation 



In recent years, the level-set framework of 21 has gained wide popularity in 
the area of image segmentation for a number of its remarkable advantages. 
In particular, level-set methods allow one to perform numerical computations 
involving curves and surfaces without the need to explicitly parameterize these 
objects. Thus, in the current setting, the subset f2i„ can be defined as 

0„ = {x e 17 I (^(x) < 0} , (5) 

where : — > M is a level-set Junction^ whose zero level-set {x e J7 | (/3(x) — 0} 
is used to implicitly define the related active contour. 
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Expressing flin in terms of its associated ip leads to a different definition of 
the Bhattacliaryya coefficient, which now becomes a function of (/?, viz. 



N N 



) = Yl ^fc(v) "II/ VP-mizk)pizk I -f) dzk, (6) 

fc=l k=l-' 

where the densities p{zk \ Lp) can be computed according to 

, , , Jfc(x))H(-^(x))rfx 
A2^(-</'(x))rfx 

with T-Llx) = {x)+ standing for the Heaviside function. Subsequently, the prob- 
lem of finding an optimal flin can be supplanted by an equivalent problem of 
finding an optimal level-set function if* as given by 

Lp* = a,Tg max B{ip), (8) 

where $ denotes a set of functions to which ip can be formally ascribed. Thus, 
for example, it is common to define $ to be a set of signed distance functions 
— the choice which leads to particularly efficient numerical implementation of 
(H [21] 



Due to the absence of a closed form solution to the above problem, a nu- 
merical scheme for maximizing ([6]) is needed. In the case at hand, a standard 
approach to the solution of Q is by means of a steepest ascent procedure which 
prescribes approximating ip* as a stationary point of the sequence of solutions 
produced by the gradient flow 

dt ~ ' ^ ' 

where a virtual (iteration) "time" t is introduced to allow the level-set function 
to evolve into a family (^(x, t), and 6B{ip)/5(p stands for the first variation of B 
w.r.t. the level-set function. Particularly, straightforward computation leads to 

^ = ^.||V^II, (10) 
with ||V(p|| denoting the (Euclidean) norm of the gradient Vtp of ip, and 

^b(x)==— ^afc Bfe(v3)- r{zk\p)*Kizk) 



2A- 

k= 



(11) 



where A := n{-(p{x))dx, r{zk \ ip) := y/pm{zk)/p{zk \ ip), at := Ilili.z^^fc Bi{ip), 
and the asterisk denotes the operation of linear convolution, which can be per- 
formed using any fast algorithm. 

It should be noted that the estimation of p* has been so far performed based 
on regional features of the observed image /. This approach is particularly use- 
ful in the cases when the objects of interest do not possess well-identifiable 
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edges 12 . However, in situations when the edges can be detected, it would be 



an omission not to consider such an important source of information. Techni- 
cally, edge-related information can be incorporated via using the framework of 
geodesic active contours 
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31 . In this case, the optimal level-set function (p* 



is estimated according to 



if* (x) = arg max 



[ g(x)||VH(<^(x))||dx 
Jn 



(12) 



where a > is a regularizing constant]^ and g : il — > M+ is an edge detector 
function, which can be defined as, e.g., g = {1 + ||V/||^)^^ with / being a 
smoothed version of /. It is straightforward to show that the gradient flow 



corresponding to (12 1 is given by 



dt 



= [ aVs + div g 



(13) 



It is worth noting that the Vb term in (131 attempts to find a region r2,;„ in the 



image domain whose empirical density P(z | Qin) closely matches the model 
density Pm{z) in terms of the Bhattacharyya metric. The second term, on the 
other hand, attempts to match the boundary of fli^ to the strong edges of /, 
while maintaining minimal curvature of the related active contour. 



2.3 Tracking of shape features 

The optimization ( |T2| ) is particularly useful in the cases when the edges of the 
object(s) of interest can be well-defined. Unfortunately, due to the possibility 
of background clutter and the presence of strong noise, using the edge infor- 
mation alone may not result in a sufficient regularization force, in which case 
the gradient flow of ( [l3| needs to be subjected to additional constraints. A 
common solution to the aforementioned problem lies in the concept of shape 
priors [3|[l4|[T5j[T8j[l9] . In this paper, as a means to distinguish the proposed 
"weak" shape priors and the above mentioned methods, the latter are consid- 
ered as based on "strong" shape priors, which impose a substantial restriction on 
the evolution of active contours. As will be demonstrated shortly, while highly 
desirable in certain scenarios, such "strong" priors can bias the segmentation 
towards a suboptimal (or even useless) solution in situations when the training 
set consists of a relatively small number of examples. The "weak" shape pri- 
ors, on the other hand, are less susceptible to the above limitation which, in 
combination with their property of being invariant under a group of geometric 
transformations, makes the proposed priors a viable alternative to the existing 
shape models. 

As a means to create a "weak" shape model based on the distribution track- 
ing framework, one must first extract a morphological feature (or a set thereof) 
which will be used to characterize the boundary of the object of interest. To 

^In the current paper, the value of a has been set to be equal to 2. 
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this end, let ip"^ be the level-set function corresponding to a known n^, and 
let Tm = {x g I (^'"(x) = 0} be the boundary of ^Im- In the current paper, 
the "weak" shape priors are constructed using the curvature of Tm as a mor- 
phological feature. Note that the curvature has been chosen primarily for its 
property of being invariant under the group of Euclidean transformations - the 
property that allows the algorithm to forgo the preprocessing step of alignment 
and registration via rotations and translations. It should be noted, however, 
that the proposed approach is by no means restricted to the tracking of curva- 
ture alone, as the very same methodology could be applied to other geometric 
descriptors of curves, should one require invariance under a different class of 



the 



geometric transformations 32 33 . 

In the case when (f"^ is defined to be a signed distance function 
values of the curvature of r„i can be obtained from the values of the curvature 
Km of the level sets of (p™ computed according to 



-div 



(14) 



Thus, given a training set of segmented images, the model pdf of the object 
curvature can be estimated as given by 



Cm[ri)=£ 



J^S,i^"'iic))K{7j-Kmi^))dy 
^<5,(^™(x))dx 



(15) 



where £ denotes the operator of assemble average as before, and denotes a 



smoothed delta function, which can be defined, e.g., as 12 



(2e)-i(l-hcos(7ra;/e)), < e 
0, otherwise. 



(16) 



Note that the role of 6^ in (151 is to "weight" the domain so as to make 
the estimated pdf depend on the values of Km in close proximity of the active 
contour. (In our experimental study, the parameter e is equal to 2.) It is 



also worthwhile noting that the estimate in (151 is conceptually analogous to 
weighted kernel density estimation [26| . Just as in the case of the "photometric" 
densities, the model pdf Cmiv) is supposed to be pre-computed and stored 
before the actual segmentation is initialized. 

Given a level set function ip, the empirical probability density C{ri \ p) of 
the curvature of its associated active contour can be estimated according to 



Civ I V) 



S^ipj^)) K{v + div{Vy.(x)/|| Vy(x)||}) dy 
^(5,(V3(x))dx 



(17) 



Subsequently, the similarity between the model and empirical densities can again 
be assessed using the Bhattacharyya coefficient, which in this case is given by 



BAp) - / VCrn{v)C{v\p)dr,, 



(18) 
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where subscript c has been added as a means to discriminate between the above 
coefficient and the one in 

Analogously to the tracking of photometric variables, we propose to track 
the empirical density C{r] \ (f) by forcing it to be as similar as possible to the 



model density Cm{ii) through maximization of Bdip) in ( 18 ). The maximization 



can be achieved using the same tool of steepest ascent in the direction of the 
first variation of Bc{^) computed w.r.t. ip. Under the assumption on if to be 
a signed distance function (which implies ||V</5(x)|l — l,Vx), the above first 
variation is given by 



K(x) 



5B,{ip) ^1 
2 



Sip 



(19) 



where the prime stands for differentiation w.r.t. 77, k(x) = — div{V(p(x)/||V(^(x)||}, 
L{ri) = y/ Cm(?/)/C(?7 I if), and A is the operator of Laplacian. Subsequently, 
the optimization problem of ( 12 ) can be extended to incorporate the additional 
shape information, which results in a new optimal ip* computed according to 



arg max 



aB{ip)+(3B,{ip)- / .9(x)||VH((^(x))||rfx 



(20) 



where $ denotes the set of signed distance functions (which can be identi- 
fied with the set of solutions of the eikonal equation ||V(^|| — 1, subject to 
</3(x)|^^g^ — 0), and /3 > is another regularization parameter, which is set 
to be equafto 5 in the experimental study of this paper. Finally, the gradient 



flow that corresponds to the maximization problem ( 20 1 is given by 



dt 



V(/3|| a Vb + div g 



+ l3Vc. 



(21) 



In the present study, the gradient flow in (21) has been implemented using an 



implicit discretization scheme based on the method of additive operator splitting 
(AOS) as detailed in [ij Ch. 3-4]. 



3 Numerical Considerations 
3.1 Regularization of curvature 

In numerical computations, the values of the curvature k of the level-sets of ip 
are commonly computed according to 

( Vip \ 2ip^iPyip^y- ip^^ifil- iPyyipl 
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Before the regularizalion After the regularization 




Figure 1: (Subplot A) A level-set function and its corresponding zero level-set 
before regularization. (Subplot Al) A zoomed section of the level-set function in 
Subplot A. (Subplot A2) The curvature of the level-sets in Subplot Al. (Subplot 
B) The regularized version of the level-set function in Subplot A. (Subplot Bl) 
A zoomed portion of the regularized level-set function in Subplot B. (Subplot 
B2) The curvature of the level-sets in Subplot Bl. 
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where subscripts x and y denote partial differentiation along the respective 



directions. Since the partial derivatives in ( 22 ) are standardly approximated by 



means of finite-support numerical schemes, the practical values of k should be 



expected to be flawed by quantization noises 34 . Consequently, the empirical 
pdf of K will be a poor estimate of the original curvature, if the effect of the noises 
is not appropriately dealt with. The problematic aspect of the above argument 
is exemplified in Fig. [l] where Subplot A shows a signed distance function 
(visualized as a grayscale image) corresponding to a circle that is depicted in 
yellow in the same subplot. In addition, Subplot Al shows a zoomed section 
of Subplot A (where the effect of quantization can be clearly observed), while 
Subplot A2 superimposes the circle over an image of its related curvature k. 
One can see that the approximated values of k suffer from spurious variations, 
which contradict the theoretical behavior of k. 

In this paper, to alleviate the problem of irregularity of the discrete ap- 
proximation of K, we propose to subject the level-set function f{-,t) to a pre- 
smoothing procedure, before f{-,t) is used to estimate the pdf of k according 



to (17). Specifically, the pre-smoothing can be achieved via anisotropically dif- 

ig equation [35] 

= div(I?(x)V¥>(x,T)), (23) 



fusing (p by means of the following equation [35| 
dfix, t) 



dr 

where -D(x) S M^^^ denotes a diffusivity tensor at x, and t denotes an artificial 



diffusion time (which is not to be confused with the one in (21)) 



The smoothing effect of the diffusion ( 23 ) is controlled through the spectral 
properties of -D(x) which, in the case at hand, should be defined in such a way 
that the smoothing will propagate along the directions tangent to the level-sets 
of (p. Note that diffusing in the tangential direction allows preserving the shape 
of the active contour, while effectively suppressing the spurious irregularities in 
if which have been induced by discretization. In particular, the above objective 
can be achieved via setting D(x) = 7t/i(x)?/i(x)-^ + -y2(x)w2(x)^, with 7^1, 
and with wi(x) and W2(x) being two unit vectors pointing in the directions 
parallel and perpendicular to V(/3(x), respectively. Namely, 'iTi(x) || V(^(x) and 
1/2 (x) _L V(^(x) at every x g fi. 



In the experiments reported in this paper, the diffusion in (23) was per 



formed using the AOS scheme as detailed in 35 , with 7 = 0.01, and with the 
number of iterations and the diffusion time-step being equal to 4 and At = 5, 
respectively. The advantage of the proposed preprocessing stage can be appre- 
ciated via observing Subplots Bl and B2 of Fig. [l] which show a regularized 
version of the level-set function of Subplot A and its associated curvature, cor- 
respondingly. One can see that the preprocessing brings the discrete values of 
K to a closer correspondence with its theoretically predicted values. 

3.2 Numerical algorithm 

For the sake of reproducibility, a pseudocode of the proposed algorithm is sum- 
marized in Algorithm 1 below. Input into the algorithm consists of the original 
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image I, a model distribution of its photometric features Pm, a model curvature 
distribution Cm, as well as of an initial level-set function ipQ. The output of 
the algorithm consists of an optimal level-set function ip*, whose zero level-set 
defines the boundary of the object of interest. 



Algorithm 1 Proposed segmentation procedure 



0) 



given: /, P„, C,„, (po{-) = ip{-,t 
preset: Ai = 5, a = 2, /3 = 5, t = 
compute: {Jk}k=i and an edge detector function g 
vifhile S > 10~^ do 

Diffuse (/j'^*-' using (23 1 to result in i^^*^ 



Compute K = -div{V(^(*Vl|V<^^*^|!} using (22 1 



N 



rr 



Compute {p{zk \ using (|7) and C(?7 

Compute Vb using and using ( 19 1 



,(t) 



^(t+i) ^ ^{t) + At (a + /3 Fc) 
^(t+i) ^ AOS ((p(*+i),g,Ai) 
Redistance by fast marching 



) using ^ 



(5-4= 1193 

end while 
return Lp* 



(t-i-i) 



(*)| 
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4 Experimental Results 

4.1 Segmentation evaluation metrics 

In this section, to qualitatively compare the performance of the proposed and 
reference segmentation methods, a number of comparison metrics are employed. 
Specifically, given the true boundary of a desired object and its estimate, both 
contours can be sampled at the points defined by their intersection with L 
rays emanating from a common barycenter of the contours (see Subplot A of 
Fig.[2]for an example). Subsequently, along the directions of the rays (defined by 
angles 9i = iTxijL, with i = 0, 1, . . . , L), the signed distances {d(Qi)\^^^ between 
the original and optimal points can be computed and used to compute mean 
difference (MD), mean absolute difference (MAD), and maximum difference 



(MAXD) defined as 36 



MD = 



\Y,d{B.,\ (24a) 

1 ^ 

MAD= -^|d(^,)l, (24b) 



i=l 



MAXD= max {|(i(6'i)|}- (24c) 

l<z<L 
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Figure 2: (Subplot A) Local difference d{9i) between the original and opti- 
mal boundaries; (Subplot B) Different regions within the original and optimal 
segmentations. 



Whereas MD provides useful information on the size and direction of the seg- 
mentation bias, MAD and MAXD quantify the overall segmentation error in 
terms of the average and maximum value of {\d{6i)\}f^i, respectively. 

Additional performance measures can be derived based on area-based met- 
rics, which (as opposed to MD, MAD, and MAXD) are insensitive to the local 
geometry of segmentation boundaries. Such metrics can be defined using the 
nomenclature of detection theory, in which case they are computed based on the 
differences between the original and optimal segmentations viewed as subsets 
of the image domain D,. Particularly, let true positive^ false positive, and false 
negative subsets of fl be defined as shown in Subplot B of Fig. [2] with their cor- 



responding areas denoted by TP, FP, and FN, respectively. Then, following 36 
one can define three area-based metrics as given by 

TP 

SEN= , (25a) 

TP -I- FN' ^ ^ 

FP -I- FN 

*° = TP + FP + FN ' <^^'=' 

With TP -I- FN being equal to the area of the original segmentation, the SEN, 
ACC, and AO measures quantify the sensitivity, accuracy and area-overlap of 



the optimal segmentation, respectively 36 . It is worthwhile noting that, in the 



case of errorless segmentation, all the measures in ( 25 ) attain their maximum 
value of 1. Thus, in general, the higher the values of SEN, ACC, and AO, the 
better the performance of an image segmentation algorithm under consideration. 
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Finally, the segmentation accuracy can be also assessed in terms of the the 
standard mean squared error (MSE) criterion. Specifically, let F be the M x N 
array of the true (binary) segmentation mask, and F be its estimate. Then, the 
MSE is defined as 



MSE-^ ^ 5^ |nx,y)-F(x,y)|'. (26) 

x=0 y=0 

4.2 Experiments with MPEG shape data 

In this section, shape data from the MPEG-7 core experiments CE-Shape-1 
database ^37j is exploited to validate the performance of the proposed seg- 
mentation method. In particular, the "face" and "teddy" data-sets are used 
throughout this section. From the set of 20 images in each data-set, a total 
of 20 leave-one-out validation experiments were performed for which 19 images 
were used for training, while the remaining image was used for testing. Hence, 
the comparative figures summarized in Table [l] below have been obtained via 
averaging the results of 20 independent trials. 

Each segmentation experiment was preceded by a training stage, during 
which a total of 19 original images in each data-set were used to learn the 
model densities Pm(z) and Cm(?7), while the remaining 20"^ image (used for the 
actual segmentation) was intentionally left out[^ The photometric pdf Pm(z) 
was computed based on two image features, viz. the values of a data image / 
as well as of its linearly smoothed version J, i.e. Ji — I and J2 = /. The data 
images, in turn, were obtained from the original images by contaminating them 
with both horizontal and vertical line clutter, supplemented by the addition of 
white Gaussian noise to further compromise the edges (see Subplots C in Fig. |3] 
and Fig. |4]for some typical examples of resulting images). 

As shown in Fig. [3] and Fig. |4] the horizontal and vertical line clutters along 
with Gaussian noise contamination result in missing edges in both the "face" 
and "teddy" images. Needless to say, these missing edges would cause problems 
for segmentation algorithms which do not rely on prior shape information. The 
proposed algorithm, on the other hand, does not have any difficulties converging 
to a close approximation of the true shape. The observed regularization occurs 



due to the Vc force as defined by ( 19 ); if the active contour were to "leak" out at 
the cluttered areas, the empirical pdf C{r] \ ip) would increasingly deviate from 
the model distribution Cmiv)- Thus, even though the shape prior "encoded" 
in terms of a pdf distribution is considered a "weak" prior, it appears to be 
sufficiently restrictive to effectively regularize the segmentation. 

To prove that the proposed method constitutes a viable and useful alterna- 
tive to existing segmentation techniques, its performance should be compared 
with that of a well-established reference algorithm. Selecting such an algorithm 
is obviously a difficult task in view of the availability of a plethora of powerful 



■^Note that, in this case, the expectations in |TJ and ^ are approximated by the assemble 
average of 19 probabiUty densities. 
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Figure 3: Segmentation of a "face" image with compromised edges. (Subplot 
A) Original binary image; (Subplot B) Original binary image with added clutter; 
(Subplot C) Cluttered image contaminated with Gaussian noise; (Subplot D) 
Initial segmentation; (Subplot E) Final segmentation. 



Table 1: Segmentation results pertaining to Section 4.2 





Face 




Teddy 


Method 


Weak Priors (icr) 


Strong Priors (±cr) 


Weak Priors (± 


cr) Strong Priors (±cr) 


MD 


0.535 ±0.189 


0.548 ± 0.823 


0.758 ±0.568 


-1.45 ± 1.06 


MAD 


2.240 ±0.101 


2.410 ± 1.230 


2.970 ±0.383 


4.02 ± 1.18 


MAXD 


4.41 ± 0.39 


3.91 ±0.44 


4.66 ±0.35 


4.85 ±0.25 


SEN 


0.974 ±0.001 


0.971 ±0.016 


0.942 ± 0.005 


0.97 ±0.015 


ACC 


0.959 ±0.002 


0.953 ±0.022 


0.936 ±0.005 


0.93 ±0.024 


AO 


0.960 ±0.001 


0.954 ±0.021 


0.937 ±0.005 


0.93 ±0.022 


MSE 


0.017 ±0.001 


0.021 ±0.009 


0.020 ±0.002 


0.023 ± 0.008 
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Figure 4: Segmentation of a "teddy" image with compromised edges. (Subplot 
A) Original binary image; (Subplot B) Original binary image with added clutter; 
(Subplot C) Cluttered image contaminated with Gaussian noise; (Subplot D) 
Initial segmentation; (Subplot E) Final segmentation. 
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Figure 5: (Subplot A) Segmentation using "weak" shape priors; (Subplot Al 
and Subplot A2) Zoomed-in sections of the segmented "face" image in Subplot 
A; (Subplot B) Segmentation using "strong" shape priors; (Subplot Bl and 
Subplot B2) Zoomed-in sections of the segmented "face" image in Subplot B. 



candidates. However, as the primary contribution of the present paper is in 
introducing a new type of shape priors, it seems reasonable to require the refer- 
ence method to be identical to the proposed one, except for the shape modeling 
part. Consequently, in this paper, the method of fl9' was used for comparisons. 



Note that the only difference between the proposed method and the one in 19 



is that the shape modeling employed by the latter is based on the PCA-based 
approach of 14 , which will be referred below to as using "strong" shape priors. 

Analogously to the above case, the prior models of the reference method 
were learned using the "leave-one-out" approach, with 19 images used for train- 
ing and one image used for validation. A typical segmentation produced by the 
proposed method is depicted in Subplot A of Fig. [Sj while Subplot B of the 
same figure shows the segmentation produced by the reference method. One 
can see that the proposed method is able to accurately delineate the "face" 
image as can be evidenced in the zoomed-in portions of the image depicted in 
Subplot Al and A2. On the other hand, Subplots Bl and B2 reveal estimation 
bias in the results obtained using the reference method. Similar segmentation 
results for the "teddy" image are depicted in Fig. |6] Here again the segmenta- 
tion obtained based on the "weak" shape model outperforms the segmentation 
obtained using the "strong" priors. The above observations and conclusions are 
further supported by the quantitative comparisons of Table [l] Particularly, in 
the set of "face" test images, the "weak" priors outperform the "strong "priors 
in all the performance metrics except for the MAXD measure. Similarly, in the 
set of "teddy" test images, the "weak" priors outperform the "strong" priors in 
all the performance metrics except in the SEN category. 
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Figure 6: (Subplot A) Segmentation using "weak" shape priors; (Subplot Al and 
Subplot A2) Zoomed-in sections of the segmented "teddy" image in Subplot A; 
(Subplot B) Segmentation using "strong" shape priors; (Subplot Bl and Subplot 
B2) Zoomed-in sections of the segmented "teddy" image in Subplot B. 
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4.3 Sensitivity to shape outliers 

The effectiveness of PCA-based prior modeling [3 14 19 is known to be depen- 
dent on tlie size of training sets used, with larger sets resulting in more robust 
and reliable segmentation. However, practical situations are common when ac- 
cess to an abundance of training shapes may not be available. Consequently, 
in such cases, the performance of PCA-based shape modeling is likely to dete- 
riorate. Furthermore, as a general rule, the smaller the size of a training set is, 
the poorer is the robustness of PCA-based modeling towards erroneous training 
samples (outliers), which should be expected in practice. The proposed "weak" 
modeling, on the other hand, are substantially less sensitive to the influence of 
outliers, and it requires substantially fewer training shapes to produce reliable 
shape models as demonstrate by our next experiment. 

In order to quantitatively compare the robustness of the "weak" and "strong" 
shape models, the same data-sets consisting of "face" and "teddy" images was 
used. In this case, however, the size of the training sets was reduced from 19 to 
4 sample images. Moreover, each training set was "spoiled" by supplementing 
it with an outlier image of an elliptic blob (whose area was set to be approxi- 
mately equal to the area of a "face" ) . Subsequently, a total of 5 training images 
were used to learn the "weak" and "strong" shape models of the test images. 
Note that, analogous to the case of Section |4[2j the data images used in actual 
segmentation have been corrupted by both vertical and horizontal linear clutter 
and white Gaussian noise. 

Fig. [T] demonstrates a typical result of the current experiment. Particularly, 
Subplots Al and A2 of the figure show the segmentations computed using the 
proposed "weak" priors, in which case the presence of outliers does not seem to 
have any effect on the performance of the algorithm. The PCA-based "strong" 
modeling, on the other hand, results in useless segmentation (as depicted by 
Subplots Bl and B2 of Fig. [?]) which supports the above concern regarding its 
being susceptible to erroneous training samples. This fact is also supported by 
the quantitative metrics of Table [2] which clearly demonstrate the superiority 
of the proposed method over the reference one. 



Table 2: Segmentation results pertaining to Section 4.3 



Method 


Face 

Weak Priors (icr) Strong Priors (±cr) 


Weak Priors (± 


Teddy 

cr) Strong Priors (±cr) 


MD 


0.571 ±0.178 


1.46 ±0.872 


0.784 ±0.556 


-12.4 ± 1.13 


MAD 


2.25 ±0.110 


6.15 ±0.965 


2.99 ±0.378 


18.5 ± 1.01 


MAXD 


4.53 ±0.39 


7.45 ±0.72 


4.66 ±0.342 


5.42 ±0.68 


SEN 


0.974 ±0.001 


0.916 ±0.017 


0.941 ±0.005 


0.91 ±0.012 


ACC 


0.958 ±0.002 


0.884 ±0.018 


0.935 ±0.005 


0.725 ±0.019 


AO 


0.960 ±0.003 


0.888 ±0.017 


0.936 ±0.004 


0.768 ±0.014 


MSB 


0.017 ±0.001 


0.072 ±0.099 


0.021 ±0.002 


0.087 ±0.006 
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Figure 7: (Subplot Al and Subplot A2) Segmentations results for "face" and 
"teddy" images obtained with the proposed algorithm. (Subplot Bl and Sub- 
plot B2) Segmentation results for the same images obtained using the reference 
method of flOl . In both cases the training data-sets consist of four correct images 
and one outlier. 
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4.4 Partial occlusions 



The experiments described in the preceding sections of the paper have been 
based on synthetic data images. In this section, the performance of both the 
proposed and reference algorithms is assessed using the real-hfe, transrectal 
ultrasound (TRUS) images. In general, ultrasound imaging is known to be an 
integrate part of modern image-based diagnostics, where it is distinguished from 
alternative imaging modalities because of its high benefit-to-cost ratio. Unfor- 
tunately, numerous advantages of ultrasound imaging (e.g. non-invasivenss, 
portability, cost efficiency, etc.) are counterbalanced by its relatively poor reso- 
lution and contrast, as compared with MRI and X-ray CT. Moreover, ultrasound 
images are also known to suffer from refraction and shadowing artifacts, which 
effectively occlude the boundaries of studied organs, thereby substantially com- 
plicating the process of their automatic segmentation. These artifacts can be 
observed in Subplots A1-A2 of Fig. [8] which show three examples of TRUS im- 
ages, while Subplots B1-B3 of the figure show the shapes of the corresponding 
prostate glands which have been manually delineated by a radiologist at the Uni- 
versity Hospital of London (Ontario, Canada). In what follows, the results of 
manual segmentation will be used as a reference, against which the performance 
of segmentation methods will be compared. 

In the experiments of this subsection, a set consisting of six TRUS images 
was used for training and validation. The photometric features used by the 
proposed and reference algorithms were defined to be the gray-level values of the 
TRUS images together with the values of their de-speckled versions computed 
by the SRAD filter of 38 . As in the previous experiments, curvature has been 
employed once again as a geometric shape descriptor. For the images in Subplots 
A1-A3 of Fig. [8j their corresponding segmentation results obtained using the 
proposed and reference algorithms are shown in Subplots C1-C3 and Subplots 
D1-D3, respectively. Despite the apparent similarity of these segmentations, 
a closer inspection reveals that the proposed algorithm provides results which 
better match the expert delineation as compared with the results computed by 
the reference method. This observation is further supported by the figures in 
Table [3l which compares the algorithms in terms of the metrics of Section |4.1 



Table 3: Segmentation results pertaining to Section |4.4| 



Method 


Weak Priors ±cr 


Strong Priors ±cr 


MD 




4.60 ± 1.96 


1.55 ±2.70 


MAD 




5.39 ± 1.51 


6.53 ±2.10 


MAXD 




7.12 ±0.708 


9.87 ±0.849 


SEN 




0.939 ±0.032 


0.902 ± 0.039 


ACC 




0.873 ±0.029 


0.850 ±0.052 


AO 




0.881 ±0.027 


0.858 ±0.047 


MSB 




0.031 ±0.008 


0.0527 ±0.019 
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Figure 8: (Subplots A1-A3) Original TRUS images of the prostate; (Subplots 
B1-B3) Manual segmentations of the prostates; (Subplots C1-C3) Segmenta- 
tion results produced based on the "strong" shape model; (Subplots D1-D3) 
Segmentation results produced using the proposed approach with /3 — 2. 
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5 Discussion and Conclusion 



The image segmentation algorithm proposed in this paper can be seen as an 
extension of the distribution tracking approach of 23 , similarly to which the 
proposed method performs segmentation via minimizing a distance between the 
empirical and model pdf's of image features. Learning the model densities is 
based on training examples, which consist of a set of representative images and 
their corresponding segmentation masks. In addition to their standard use as 
identifiers of the object class, the masks can be also used to estimate the pdf 
of geometric parameters of the object boundaries. This allows the approach 
of 23 to be extended to tracking the "morphological distributions" , in which 
case the optimal active contour is required to minimize a distance between the 
model pdf of the geometric parameter(s) and its empirical counterpart. In the 
current paper, the Bhattacharyya coefficient has been employed to assess the 
distances between the model and empirical distributions for both photometric 
and morphological features. Needless to say, the chosen distance can be re- 
placed by other metrics (such as the KuUback-Leibler divergence or the Fisher 
discriminant [30]), should there be a rationale behind such a modification. 

It goes without saying that, given a pdf of a geometric parameter (such as, 
e.g., curvature), there exists a myriad of possible shapes whose boundaries will 
have their empirical pdf's identical to the one above. For this reason, the pro- 
posed approach to shape modeling has been coined as "weak" - the term which 
is meant to accentuate the minimally restrictive nature of corresponding shape 
priors. It has been proven via a series of numerical experiments that the "weak" 
modeling is nevertheless capable of effectively regularizing the convergence of 
active contour under conditions of cluttered and occluded object boundaries. 
Moreover, as opposed to the PCA-based modeling 3 p4p!9 , learning the "weak" 
priors can be performed using relatively small training sets without noticeably 
compromising the performance of resulting segmentation. This makes the pro- 
posed method particularly useful in situations when one does not have access 
to an abundance of training shapes. 

In practice, the images comprising a given training set are normally reg- 
istered to a common position, orientation, and scale before a specific shape 
model can be learned. Within an actual image, on the other hand, the object 
to be segmented can, in general, appear to be unaligned with respect to such 
model. For this reason, many segmentation algorithms are bound to use an 
image registration procedure as a means to bring an evolving active contour 
into correspondence with its model before the discrepancy between them can be 
assessed. Needless to say, the registration can substantially increase the overall 
cost of image segmentation, which is clearly a disadvantage. The "weak" shape 
modeling, on the other hand, is free of the above limitation. Indeed, by choosing 
the shape descriptor to be invariant under a group of geometric transformations, 
one can effectively forgo the image alignment stage (as long as the latter is based 
on the assumed type of transformation). Note that, in the present work, curva- 
ture has been used as a shape parameter in the "weak" modeling, which makes 
the resulting segmentation invariant under the group of Euclidean transforma- 
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tions. Obviously, the curvature would have to be replaced by a different shape 
descriptor, should one require the invariance under a different (larger) class of 
geometric transformations 32 33 . 
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